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SUMMARY 

This Interim report consists of two reports: "Space Radiation Effects 
on Si APDs for GLAS" and "Computer Simulation of Avalanche 
Photodiode and Preamplifier Output for Laser Altimeters." The former 
contains a detailed description of our proton radiation test of Si 
APDs performed at the Brookhaven National Laboratory. The latter 
documents the computer program subroutines which we wrote for the 
upgrade of NASA's GLAS simulator. 
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Abstract 

The test results are reported of proton radiation damage of Si APDs for 
use in the GLAS laser altimeter and the cloud lidar. The APD bulk leakage 
current increased at 0.26 fA/rad, or about 1600 dark counts/s per rad at -10°C 
under 16.2 MeV protons. The space radiation damage is expected to cause the 
altimeter receiver noise spectral density to double after 3 krad at the end of 
the mission. The radiation damage to the Si APD single photon counter for 
the GLAS cloud lidar is expected to be much more severe and the receiver 
performance will be seriously affected, especially during nighttime operation. 
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1 Introduction 


The Geoscience Laser Altimeter System (GLAS) consists of a laser altimeter at 1064 
nm wavelength and a cloud lidar at 532 nm wavelength. A Si avalanche photodiode 
(APD) will be used as the altimeter photodetector. A photomultiplier tube or a 
Geiger mode Si APD single photon counter is currently being considered as the cloud 
lidar photodetector. It is estimated that GLAS will accumulate about 3 to 4 krad 
of space radiation over the 5 year mission with 7 mm thick aluminum shielding [1]. 
The type of radiations are primarily protons in the South Atlantic Anomaly and solar 
protons at the poles of the earth geomagnitosphere [2]. The energy spectrum of the 
protons is centered at about 50 MeV and extends from 1 MeV to several hundred 
MeV. We recently measured proton radiation damages of Si APDs for the GLAS 
mission. The results shows that Si APDs in linear mode for the GLAS altimeter will 
suffer significant but tolerable radiation damage after 3 krad total dose. However, 
Si APDs in Geiger mode as single photon counters are shown to suffer significant 
damages after as little as 100 rad and started to saturate after 3 krad at — 10°C, and 
200 rad at room temperature. 

Due to their mass and charge, protons cause both ionization damages to the sur- 
face insulator layer and displacement damages to the crystal lattice in semiconductor 
devices [3]. Ionization damage is the major concern for surface devices, such as MOS- 
FETs and CCDs. Displacement damage effectively adds defects to the semiconductor 
material and causes degradations in the device characteristics which are critically de- 
pendent on the number of defects in the material. Ionization damage is often tested 
with gamma rays because they are easy to obtain and has the same effects with re- 
spect to ionization damage. Displacement damage is often tested with neutrons which 
mainly cause displacement damage. Raymond and et al. showed that 1 rad of 1-100 
MeV protons are equivalent to a fluence of 1 to 1.5 x 10 ‘ j cm? 1 MeV neutron with 
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respect to displacement damage [4]. 

Previous studies of radiation damages of PIN photodiodes and APDs used either 
gamma rays or electrons ( 60 Co) or neutrons as the radiation source [5]-[13], Based 
on those studies, it is estimated that the Si APDs for the GLAS laser altimeter will 
suffer a considerable increase in the bulk leakage current after 3 krad though the 
resultant degradation in the overall receiver performance is still tolerable. On the 
other hand, the damage to the Geiger mode APD single photon counter for the cloud 
lidar receiver is expected to become catastrophic after several hundred rad, because 
each electron of the bulk leakage current will be counted as a dark count. Other 
APD parameters such as the quantum efficiency, the gain, the device speed, and the 
breakdown voltage are not expected to change significantly for the relatively low dose 
that the GLAS receiver is subject to [8] [9]. 

There has been no direct measurement of proton radiation damage of Si APDs. 
It is of interest to compare the measurement results with the estimates based on 
the neutron radiation damage test. It is also of interest to test the proton-neutron 
equivalence relationship by Raymond for photodiodes, especially APDs which are 
normally operated at a much higher electrical field intensity than other types of 
semiconductor devices. 

The rest of this report describe the test setup, the measurement procedure, the 
test results, and the estimated space radiation damage for the GLAS laser altimeter 
and cloud lidar receivers. 

2 Experiment Setup 

The proton radiation source used was the Tandem Van de Graaff accelerator at the 
Brookhaven National Laboratory [14]. The tests were conducted on April 25-26, 1994. 
Three proton energy levels were chosen: 14.8 MeV, 22 MeV, and 28 MeV. The flux 
was set to 10 6 — 10' protons/cm 2 sec. 
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The test samples consisted of nine EG&G C30902STC devices, three for each 
proton energies. These devices are used by EG&G in their single photon counting 
module SPCM-100-PQ. Each device contained a C30902S Si APDs on a small thermal 
electric cooler, all in a TO type package. The diameter of the APD active area was 
0.5 mm and the thickness of the active volume was 27 ^m. The active volume was, 
therefore, 5.3 x 10 -6 cm 3 . The glass window in front of the APD was 1.0 mm thick 
and was made of Bora-Silicate glass, which had a density of 2.6 and a composition of 
0(60%), Si(26%), B(9%), Na(3%), and Al(l%). The window attenuated the proton 
energy to 5.1, 16.2, and 23.4 MeV. respectively [14], The glass windows were not 
under the direct influence of the thermal electrical cooler. 

Figure 1 shows the circuit layout of the test fixture. The circuit used was a typical 
passive quenching photon counting circuit given in the EG&G data sheet [15]. The 
APDs were tested at both 22°C and — 10°C. Three APDs were mounted on the test 
fixture and were irradiated simultaneously. The proton beam had a radius of 17.5 cm 
and was uniform within 2% across the beam. The three APDs were on a 9.525 cm 
radius circle and were well within the beam. The photoelectron pulses were brought 
out via the 50D coax cable Al. Bl. and Cl. The other three coax cables, A2, B2, 
and C2, were used to carry the returns of the APD bias currents which were to 
be measured by the electrometer. This configuration helped to reduce the leakage 
current and settlement time in the low current measurements because there was no 
common mode voltage in the current measurement. As the current increased to above 
1.0 fj,A, the voltage drop across the 200AT1 resistor became significant and it should 
be subtracted from the voltage measured by the voltage meter. 

Figure 2 shows the test instrument setup. A personal computer (PC) was used 
to automatically control the measurement and to log the test data. The bias voltage 
for the APD was generated by the Analog Module 521-5 programmable high voltage 
supply and was controlled by the D-to-A output of the PC. The resolution of the bias 
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voltage output was 0.10 volts. The bias voltage was monitored independently by the 
HP3478A digital multimeter. The Keithely 617 electrometer which measure the APD 
current was put in series in the return of the APD bias supply. The voltage burden of 
the electrometer was less than 1.5 mV. The photoelectron pulses were amplified and 
then fed to the SR620 universal counter. A 50 MHz lowpass filter was used to reduce 
the noise counts caused by the circuit thermal noise. The oscilloscope displayed in 
real time the photoelectron pulses and it also measured the pulse height distribution. 
In order to avoid ground loop, the chassis of all the instruments were isolated from 
the power line ground and tied together at one point to the vacuum chamber. The 
signal ground which was carried by the shield of the coax cable Al, Bl, and Cl were 
also connected to vacuum chamber at the feed through plate. Only one APD could be 
tested at a time by connecting the appropriate type D and BNC connectors, as shown 
in Figure 2. During the irradiation, all three APDs were biased at about 120 volts 
simultaneously through the use of a special three way adapter cable (not shown). 

The dark current increase versus dose for the APDs in linear mode could also 
be determined by biasing the APD below the breakdown point. The dark current 
consisted of the bulk leakage current and the surface leakage current. The former 
passed through the AP D high field avalanche section and was the source of the dark 
counts for Geiger APDs. The latter did not pass through the avalanche region and 
might be assumed nearly unchanged as the bias voltage changed around the break- 
down point. As a worst case, one may approximate that all the radiation induced 
dark current was from the bulk leakage current. The net increase in the bulk dark 
current was equal to the difference of the measured dark currents before and after the 
irradiation divided by the average APD gain. 

Each measurement consisted two parts: determining the device break-down volt- 
age and then taking a set of measurements at a series of bias voltages around the 
break-down point. The first part of the measurement was carried out as follows. 
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(1) . The device temperature was set (manually) to either 22 °C (thermistor=11.4A'17) 
or — 10°C (thermistor=50.9AT2). 

(2) . The high voltage module output a voltage. 

(3) . The HP3478A DMM measured the bias voltage. 

(4) . The SR620 universal counter counted the number of pulses received. The count- 
ing interval was set to 0.1 seconds. The threshold level was set to 150 mV, which was 
relatively high but necessary to eliminated the noise counts due to the circuit thermal 
noise and RF pick-up from the environment. 

The process continued by repeating Steps (2) through (4) at another setting of 
the high voltage. We began at roughly 10 volts below the break-down point, and 
incremented at 1 volt steps until the number of counts exceeded 3000/s. The break- 
down voltage w ? as determined as the voltage at which the counting rate was 100/s. 
This point was obtained through linear interpolation between the two points below 
and above 100 counts/s. One could observe the pulses on the oscilloscope display 
once the bias voltage exceeded the break-down point. This definition of the break- 
down voltage might be slightly higher than the true value. Nevertheless, it served as 
a reference point in our measurements. 

Once the breakdown point was determined, the second part of the measurement 
was proceeded as follows. 

(1) . The high voltage module output a voltage. 

(2) . The HP 3478 A DMM measured the bias voltage. 

(3) . The Keithely 617 electrometer measured the current through the APD. Ten 
measurements were taken and the mean and the standard deviation were calculated. 

(4) . The SR620 universal counter counted the number of pulses received. The gate 
time and the threshold were set to 0.1 seconds and 150 mV, respectively. Ten mea- 
surements were taken and the mean and the standard deviation were calculated. 

(5) . The oscilloscope (HPo4720A) measured the mean and variance of the pulse 
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(5). The oscilloscope (HP54720A) measured the mean and variance of the pulse 
amplitude distribution. 

Steps (1) through (5) were repeated for each bias voltage setting. The voltages 
selected were 50, 150. V b -l, H-0.5, V b , K + 0.5, V b + l, V b + 2, H + 3, V b +4, V b + 5, 
V b + 6, volts with V b the measured break-down voltage. Later in the experiment, we 
also added points at V b — 10, V b — 5, V b — 4, V b — 3, and V b — 2 volts. 

The performance of those APDs as single photon counter was tested prior to the 
irradiation using a He-Ne laser (A = 632 nm). The photoelectron detection probability 
versus bias voltage was measured up to V b + 6 at 22°C. The measurement results and 
the pulse shape were close to those given in the APD data sheet [15]. The pulse width 
measured at 10% points was about 1.0 ps. The residual noise current of the entire 
setup was about 0.40 nA when the laser was off and the bias voltage was zero. 

3 Measurement Results and Discussions 

We first measured the instantaneous dark count rate of three APDs at 14.8 MeV 
energy and at a flux of about 1000 protons/ s • cm 2 . Such a flux corresponded to 
about 10 -3 rad/s after the glass window and is the highest dose rate expected for the 
GLAS receiver as it orbits through the South Atlantic Anomaly [1], There were no 
significant increase in the measured dark counts. This was expected since the number 
of protons which hit the APD active area were only 1.9 6protons/ s. The dark count 
rate of the APD. on the other hand, was typically a few hundred to a thousand per 
second at — 10°C. The radiation induced dark counts were not detectable under this 
condition unless one single proton could produce thousands of dark counts. 

The proton flux was then increased to about 10 6 /s and the measurement proce- 
dures described in the previous section were carried out with the proton beam blocked 
after a preset dose was accumulated. The first set of the three APDs were irradiated 
at 14.8 MeV and measured after 100. 200, 300, 400, 500, 1,000, 3,000, and 10,000 



rad. The second set were irradiated at 22 MeV and measured after 50, 100, 200, 
500. 1,000, 3,000, 10,000, and 30,000 rad. The third set were irradiated at 28 MeV 
and measured after 16, 31, 62, 155, 311, 933, and 3,110 rad. The measurements were 
repeated at 1 week and 1 month after the irradiation to determine the annealling 
effect. Lastly, the three APDs which were irradiated by the 22 MeV protons were 
measured one more time after being heated up in an oven at 100°C for 1 hour. 

Figures 3 and 6 show the measured dark currents and dark counts of one of the 
APD (#4, at 22 MeV) as a function of the measured bias voltage at 22°C' and -10°C, 
respectively. The other eight APDs showed similar characteristics. The decrease in 
the dark counts at 22°C as the dose increased from 3 krad 30 krad, as shown in 
Figure 4, was clearly caused by the circuit saturation when compared to Figure 6. 
As the dark counts exceeded certain level, the pulse amplitudes began to decrease 
because of a limited supply current. As a result, the number of pulses which crossed 
the discriminator threshold decreased as the actual dark counts increased. This also 
explained why the dark counts at 22°C increased from 1 week to 1 month after the 
irradiation. The slight decrease of the dark current below the break-down point at 
— 10°C (Figure 5) was believed to result from the test circuit, not the APD, because 
the measured currents were too close to the residual noise current of the setup (0.4 
nA). 

One may concluded from those figures that: (a) Significant increases in the dark 
counts began as early as 50 rad; (b) The dark counts at — 10 °C became excessively 
high, about half a million per second, after 3 krad; (c) There were little changes in 
the break-down voltages with the radiation doses; (d) The increase in the total dark 
currents below the break-down voltage was not significant until 3 krad; (e) There 
was some annealling, about 25% in the dark counts and about a factor of 2 in the 
dark current (below break-down point) for this particular device. Heating the device 
accelerated the annealling process. 
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The ratio of the dark counts at -10°C and 22°C corresponded roughly to the 
relationship given in [16], 

huik(T) <x e -°- 55eV ' kT ( 1 ) 

where T is the absolute temperature in Kelvin and k is Boltzmann’s constant. One 
may used this equation to predict the dark current and dark counts at other tem- 
peratures. As an example, the dark count rate increases by 13.7 from — 10°C' to 
22 °C. 


Figure 7 shows plots of the dark counts vs. doses using the same data set as Figure 
3. The number of the dark counts were corrected for the detector dead-time using 
the formula [17] 


1 - n(r d /T 0 ) ^ 

where h is the maximum likelihood estimate of the dark counts, n is the actually 
measured dark counts, is the dead time, and T 0 is the counting interval. The dead- 
time was taken to be r d = l.O/zs, which was roughly the same as the pulse width of the 
Geiger mode APDs in this setup. Since the dead-time correction became increasingly 
inaccurate as the value of n(r d jT 0 ) approaches unity, we eliminated in Figure 7 those 


data points with n > 0.5 x 10 6 /s. Furthermore, due to the slight variation in the 
break-down voltage at different doses, the dark counts were not all measured at 
exactly the same set of the voltage points. As a result, the voltage point set used 
in the zero dose measurement were used as the reference and the data points for the 
nonzero doses were obtained using the interpolation function of the KaleidaGraph 
software package [18]. It is shown in Figure 7 that the dark counts increased almost 
linearly with the doses and the bias voltage only affected the detection probability. 

The net increase in the dark count rate vs. proton dose may be obtained by 
first subtracting the count rate at zero rad and then dividing it by the corresponding 
photoelectron detection probability. The photoelectron detection probabilities were 
measured prior to the irradiation using a He-Ne laser and it was 18.9% for this APD 
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at 22 °C and at 5 volts above the break-down point. The measured photoelectron 
detection probabilities for lower bias voltages were not used since they were less 
accurate and depended increasingly on the threshold setting and the accuracy of 
the break-down voltage measurement. According to [15], the detection probability 
was a function of the voltage above the break-down point and independent of the 
temperature. We further assumed that the photoelectron detection probability was 
unchanged before and after the irradiation. Figure 8 plots the net dark count rate 
increase vs. dose up to 3 krad. The straight line which passes the origin was obtained 
by a least square fit. The proton radiation damage coefficient which was the slope of 
the line was 1.76 x 10 3 counts/s ■ rad. 

Table 1 shows the net dark count increase vs. dose for all the nine Si APDs. 
The 22 MeV protons caused more damage than the other two energies. The total 
detected APD dark counts at a different temperature and a different bias voltage may 
be estimated by multiplying the results in Table 1 by 


Pd(V- V b ) • exp[- 


0.55eF 1 

L V 


T{*K) 


263.5 


)] 


( 3 ) 


where Pd(V -V b ) is the photoelectron detection probability. For other models of Si 
APDs. the proton radiation damage may be calculated by scaling the above result 
for the appropriate active volume. 

We also measured the detected photon counts due to the fluorescence of the glass 
window by putting an irradiated APD face-to-face against a undamaged APD on 
the test fixture. The fluorescence photon counts were no more than 200/s measured 
immediately after 10 krad 14.8 MeV proton irradiation. Therefore, the effect of 


fluorescence from the glass window was negligible. 

It would required much higher proton doses to determine the bulk leakage current 
increase vs. dose for the APDs in linear mode when the biased voltage was about 10 
volts below the break-down point (APD gain« 100). As shown in Figures 5 and 6, 
the increase in the dark current only became apparent at voltages very close to the 
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Table 1: Test Result of the APD dark counts increase due to proton irradiation 


Proton Energy 
before Window 
after Window 

14.8 MeV 
5.1 MeV 

22 MeV 
16.2 MeV 

28 MeV 
23.4 MeV 

APD Serial # 

#1 #2 #3 

#4 #5 #6 

oo 

CO 

Acts/s ■ rad 

1540 1230 1500 

1760 2530 1600 

761 871 885 

Average 

1420 

BBEBIDi 

839 


Exclude A PD #5. 


break-down point for the proton radiation doses we used. The average APD gain at 
such a high bias voltage was too high to maintain a stable linear operation. However, 
one may assume the bulk leakage current was unchanged as the bias voltage decreased 
to below the break-down point. The bulk leakage current due to the irradiation can 
be obtained by multiplying the net APD dark count rate by the electron charge, 
e.g. 0.26 fA/rad at — 10°C and 22 MeV. The corresponding dark current at room 
temperature (22 °C) should be 3.6 fA/rad according to Eq.(l). 

The noise current spectral density due to the irradiation is given by [19] 

N rad = 2qG 2 FAI bulk (A 2 /Hz) (4) 


where q is the electron charge. G is the average APD gain, and F is the APD excess 
noise factor given by F = k e fjG + (2 — \/G){\ — k e ff ) with k e ff the APD ionization 
coefficient ratio. The corresponding noise equivalent power (NEP) is equal to \J N Tad 
divided by the APD responsivity, i.e.. 


NEP rad = 


\ A ra( i 

qG^/hf) 


\/2qG 2 F AI bulk 

qGWJTn 


{Watts fy/Hz) 


( 5 ) 


where q is the APD quantum efficiency and hf is the photon energy. 

The measurement data were several times smaller than that calculated according 
to the neutron test results and the proton-neutron equivalence relationship by Ray- 
mond. Kraner [12] showed that, for neutrons, AI bu i k = a$V with a a leakage current 
constant, $ the neutron fluence. and V the device volume. Buchinger and et al. [13] 
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showed that a = 1.7 x 10~ 16 Aj cm for this kind of APDs at room temperature. Ray- 
mond [4] showed that 1 rad protons over 1-100 MeV range had about the same effect 
in displacement damage as 1.5 x 10 7 /cm 2 neutron fluence. Therefore, the predicted 
APD bulk current due to the proton irradiation is 13.5 fA/rad at room temperature, 
which is about 4 times of what we measured (3.6 fA/rad). 

On the other hand, the radiation damage we measured was much higher than those 
obtained with gamma rays. The APD used in the Mars Observer Laser Altimeter 
(MOLA) was claimed to have no performance degradation after 100 krad total dose 
according to the specification from the McDonnell Douglas Astronautics Company. 
Conner and Guggenmos [20] also reported that, for a Si APD very similar to those 
we tested, the total dark current increased by only about 30% after 100 krad at room 
temperature. There was no further details given about those tests. 

4 Estimated Space Radiation Damage of Si APDs 
for GLAS 

4.1 The GLAS Cloud Lidar Receiver 

One of the candidate photodetector for the GLAS lidar receiver is a single photon 
counting module which contains a newly developed large area {(j> = 0.5mm) EG&G 
Shk TXI Si APD with an active quenching circuit [21]. The thickness of the device 
active region is 25 fim and the active volume is 4.9 x 10" 6 cm 3 , The active quenching 
circuit enables the photon counter to have a short dead-time, about 30 ns, which 
corresponds to a useful maximum count rate of 15 x 10 6 /s. The bias voltage is assumed 
to be 15 volts above the break-down point, which gives a photoelectron detection 
probability of 80% and a photon detection probability of 68% at A = 532 nm. For low 

dark count rates, the APD is mounted on a two stage cooler and the temperature is 
held at — 40°C. 

We assume that the proton irradiation has the same effect for the Slik APD. 
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Because of the smaller active volume and lower temperature, the dark counts of this 
Slik™ photon counter due to proton radiation is predicted to be 0.175 times those 
listed in Table 1. For example, at 14.8 MeV, the detected dark count is estimated to 
be about 250/s • rad. The dark counts of the single photon counting module will be 
about 750,000/s after 3 krad at the end of the GLAS mission, which is about 3 orders 
of magnitude higher than the dark count rate before launch. This would seriously 
affect the overall receiver performance, especially during nighttime measurement when 
the detected noise photons from moon lit clouds is only ~ 1000/s. Although the 
radiation induced dark counts are still less than the detected signal photons and 
daytime background noise photons, the effect of radiation damage has to be considered 
when designing the receiver. Further increase of the shielding thickness is not effective 
for high energy protons [1]. Further decrease of the device temperature can reduce the 
dark counts, for example, by a factor of 56 from — 40°C to — 70°C, but the detector 
subsystem may become too complicated. On the other hand, if one uses a single stage 
thermal electric cooler and the APD temperature is held at — 20°C, the dark counts 
would be 2160/s • rad, or 6.5 x 10 6 /s after 3 krad. At room temperature, the single 
photon counting module will start to saturate (15 x 10 6 /s) after as little as 200 rad. 

The transient increase in the dark count when GLAS passes the South Atlantic 
Anomaly can be neglected according to the real time measurement described at the 
beginning of the previous section. The Si APD photon counting modules are definitely 
suitable for short space missions. 

4.2 The GLAS Altimeter Receiver 

The Si APD to be used in the GLAS altimeter receiver is the same as that used 
in MOLA. Its basic structure is similar to those we tested but the active volume is 
thicker, 140 /im. Typical quantum efficiencies of 35% at A = 1.064nm is assumed, 
though 43% has been reported as the highest [22]. The ionization coefficient ratio is 
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roughly k c jj = 0.007. The diameter of the active area is <f> = 0.8mm. The active 
volume is therefore 7 x 10 -5 cm 3 . The APD is operated in linear mode at room or 
slightly higher temperature for high quantum efficiency at 1.064^tm wavelength. The 
bias voltage is set for an average gain of G ~ 200. The receiver electrical bandwidth 
is about 150 MHz. The expected increase in the APD bulk leakage current according 
to our test result is 0.041 pA/rad for 14.8 MeV protons at 22°C. 

If we assume the space radiation damage by the protons over the entire energy 
spectrum is, on the average, the same as that at 14.8 MeV, the total bulk leakage 
current due to space radiation will be 120 pA after 3 krad. The rms noise due to this 
increase in according to Eq.(4) will be 2.3pA/ Hz 1 / 2 , which is almost the same 
as the total equivalent input noise current density of the APD preamplifier module 
(2 ,6pA/ Hz 1 / 2 ) prior to irradiation. The corresponding NEP caused by the radiation 
damage is 38 fW/Hz 1//2 , or 0.47 nW over the 150 MHz bandwidth. The received 
background noise power from clouds is estimated to be about 3 nW at daytime and 
3 pW at nighttime. Therefore, the GLAS altimeter receiver will show some radiation 
related degradation, but much less severe than that of the single photon counting 
module for the cloud lidar receiver. 
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Figure 1. Circuit Layout of the Test Fixture. 
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1. Introduction 

The following is a description of the computer program which 
simulate the output signal from an avalanche photodiode (APD) and 
preamplifier in a laser altimeter receiver. The signal consists of 
random number of photoelectrons output from the APD and the noise 
from the preamplifier. The former follows the Mclntyre-Webb 
distribution [1][2][3] and the latter follows the Gaussian 
distribution. The two can be assumed independent of each other. The 
output of the simulation program consists of an array of discrete 
samples at a given sampling rate. Each sample is equal to the signal 
averaged over the sampling interval. In the calculation, we define 
the signal as the input current to the preamplifier. All the noise 
sources will be converted to equivalent current noises at the input 
to the preamplifier. 

Figure 1 shows a block diagram of the program. The input signal is 
an array of discrete sample of the received optical power. Each 
sample represents the average value of the actual signal over the 
sampling interval. The program consists of a Mclntyre-Webb random 
number generator for the APD secondary photoelectrons and a 
Gaussian random number generator for the preamplifier noise. A 
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digital lowpass filter is also included to simulate the noise filter in 
the laser altimeter receiver. 


The preamplifier noise follows the Gaussian distribution. Press and 
et al. [4] have shown a subroutine for generating standard Gaussian 
random variable. We need only to scale this standard Gaussian 
random variable, y g o, to obtain the preamplifier noise output, y g , 
according to the relationship 

y t = cry !0 +Li (1) 

where p and a are the mean and standard deviation of the 
preamplifier noise. The mean of amplifier noise can be considered to 
be zero. The variance has to be calculated in frequency domain 
because only the preamplifier noise spectral density is given. 

Since each sample represents the average signal over the sampling 
interval, i.e. 

1 r f , + Ar 

x i =—\x(t)dt (2) 

At J u 


the sampled preamplifier noise should also be the average amount of 
noise over the sampling interval. This is equivalent to have a filter 
after the preamplifier with rectangular impulse response of 
duration At and amplitude 1/Ax. The Fourier transform of this filter 
is given by 

si n_ (6 )AT /2 ) 

* coAt/2 v ' 
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It is the analog to digital sampling process that results in this 
filter. The bandwidth of the filter is roughly 1/Ax, which should be 
much wider than the actual receiver electrical bandwidth. 


The variance of the sampled preamplifier Gaussian noise can be 
calculated as 

J>,hW* (A2) (4) 


where N(co) is the one sided noise power spectrum of the 
preamplifier. 


If the preamplifier contains a FET as the front stage, the noise 
power spectrum can be written as [5] 


... . 4 KT a . _ 4KT T 4 KTT , _ x2 

N((i)) = — — — — + 2ql t + ■£— + 2— (fi)C;) 2 


R 


f 


g m R f 


Sn 


(A2/HZ) (5) 


where K is the Botzmann's constant, T a (°K) is the ambient 
temperature, Rf is the transimpedance of the preamplifier, g m is the 
transconductance of the FET, r is a factor close to unity, and Cj is 
the total capacitance at the input of the preamplifier. 


When the altimeter receiver bandwidth is below 200 MHz, one can 
often approximate the preamplifier noise as white Gaussian noise 
with the variance given as 


i 

cr = 


1 

R f 2 At 


(A 2 ) 


(6) 


where T n (°K) is the equivalent noise temperature. We used Eq. (6) to 
calculate the rms preamplifier noise in the program. 
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3 , Simulation of the Photocurrent Output from the APD 

The sampled output of an APD is a discrete random process which is 
equal to the numbers of secondary photoelectrons output by the APD 
during the sampling interval times the electron charge and then 
divided by the sampling interval. The probability of the number of 
photoelectrons has been shown to follow the so-called McIntyre 
distribution [1][2], as 

P(miX(t)Ar) = £ Py(m\n ) ■ ait)AT) " e -W f m > i 

n f 


/ > (0lA(r)Ar) = <r A( ' ,Ar 


(7) 


with Pm (min) the probability of generating m secondary 
photoelectrons in response to n primary photoelectrons, given 

f \ 


«r 


m 


+ 1 


P M (m\n) = 


1-*. 


m{m-n)\ T 


K<r m 


1 -k. 


• + 1 + n 


*!t 


by 


( 8 ) 




l + k eff (G- 1) 

n-f — = 

ra-^xG-Di 

G 


G 


where k e ff is the APD hole to electron ionization coefficient ratio, G 
is the average APD gain, and X(t) is the primary photon counting rate 
given by 

m=-i[ p ,o)+ p ,]+^ ( 9 ) 


where t\ is the APD quantum efficiency, hf is the photon energy, P s (t) 
is the input optical signal power, Pq is the background radiation 
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power assuming to be constant, lb is the APD bulk leakage current, 
and q is the electron charge. 


A very close approximation to (7) was given by Webb and et al. [3] 
forA(t)Ax»1, as 


P w (mlX(t) At) 


^Zn&FMt) At 


1 + 


1 

(m-GA(r)AT)(F- 
GFA(r) At 


a r 


exp 


(m-GA(r)Ar) 2 


2G 2 FA(r)Adl + 


(m-GA(r)Ar)(F-l) 
GFA(r)Ar 


( 10 ) 


m> A(r) A t 


where F is the APD excess noise factor which can be calculated as 


F = ^G + 



( 11 ) 


It is also known [6] that the Webb's approximation becomes 
inaccurate as A(t)Ax «1. Therefore, we choose to use Eq. (10) only 
for A(t)Ax > 1 0. 

Because the complexity of the probability distributions for the APD, 
one cannot generate the random variable as an analytical function of 
a uniform distributed random variable. Instead, the rejection method 
described in [4] has to be used. A comparison function has to be 
found when using this method. A comparison function is an upper 
bound for the desired probability distribution function and the total 
area under the curve has to be finite. The inverse function of the 
indefinite integral of the comparison function should also be known 
analytically. The area (integral) of the comparison function should 
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be as small as possible since the time required to generate a random 
number is proportional to the area. 

We simulate the APD photoelectron output for X(t)At < 10 in two 
steps. First we generate a Poisson variables to simulate the number 
of primary photo electrons. We then generate a McIntyre random 
number for each primary photoelectron according to (8). The APD 
output should be the sum of secondary photoelectrons in respond to 
all primary photoelectrons. This procedure involves calling a Poisson 
random number generator once and then calling a McIntyre random 
number generator n times with n the outcome of the Poisson random 
number generator. This method only requires a comparison function 
for the McIntyre distribution for the n=1 case. 


The Poisson random number generator can be implemented using the 
computer subroutine provided by Press and et al. [4]. The McIntyre 
random number generator can be realized similarly but with the 
McIntyre probability distribution function and an appropriate 
comparison function. We may still use a Lorentzian distribution as 
the comparison function as in [4], i.e., 


/(*) = 


(x-X :„) 2 


1 + 


( 12 ) 


but we need to find a set of values for co, xo, and ao such that 
f(m)>PM(mll) and jf(x)dx is as small as possible. As an example, we 

found co=4/G, xo=-G/50, and ao=G/3 may be used for APD gains from 
G=50 to 500. 
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The expression for McIntyre distribution (8) has to be rewritten or it 
will easily cause overflow in the computer as m become large. Using 
the property of the Gamma function, one can rewrite P(mln=1), Eq. 

(8), as 




(13) 
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For A.(t)Ax >10, one may simulate the APD output in one step using 
the Webb's approximation, Eq. (10), to save computation times. The 
comparison function is again chosen to be Lorentzian, Eq. (12), but 
with c 0 =1/G(FX(t)Ax) 1/2 , x 0 =GA(t)Ax, and ao=G(2FX(t)Ax)^2 ( W hich is 
adequate for average APD gains from 10 to 500 and average primary 
photoelectrons from 10 to 100, or a even wider range (yet to be 
tested). 

APDs also have surface leakage current which is not multiplied by 
the APD gain. This part of the APD dark current should be modeled as 

a DC current which results in an additive Gaussian noise of mean and 
variance 

n = ql, At 

a j = ql s Ar ( 1 4) 

where l s is the surface leakage current. 

4. Simulation of the Lowpass Noise FiltPr 

The lowpass noise filter is simulated using the impulse invariance 
method [7], i.e., the impulse response of the discrete time filter 
being equal to the samples of the impulse response of the 
corresponding continuous time filter, as 

h(n) = At h c (n At) ( 15 j 

where h c (t) is the impulse response of the continuous time filter to 
be simulated. This method gives the closest simulation of the actual 
filter in the receiver as long as both the signal and the filter are 
band limited and the sample rate is sufficiently high. 
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Usually, a buildable lowpass filter are specified by the locations of 
the poles and zeros each of which results from a physical capacitor 
or inductor. One can derive the system function and obtain the filter 
impulse response by taking the inverse Fourier transform of the 
filter system function. The discrete time filter impulse response 
can then be obtained according to Eq. (15). The impulse response 
function of a bandwidth limited filter is ideally infinite in time. 
However, we may truncate it as the values becomes vanishingly 
small and turn it into a finite impulse response (FIR) digital filter. 
The filter impulse response is derived separately and stored as an 
array. 

The lowpass filter output is the convolution of the input signal with 
the impulse response function, as 

y(n)= J d x(k)h(n-k + 1) (16) 

k=n-N h +1. Jr£l 

where x(n) and h(n) are the samples of the input signal and the 
impulse response of the filter. Here we assumed h(n)=0 if n<1 and 
n>Nh. The first Nh points in y(n) should be discarded because the 
input sequence is truncated, i.e. x(k)=0 if k<1. 

The choice of the filter for a given input pulse shape and noise 
characteristics may be optimized by trail and error using this 
simulation program. Usually, a Gaussian or a Bessel lowpass filter is 
used for minimum pulse shape distortion. The former is easy to 
simulate but can only be implemented approximately. The latter can 
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As an example, we implemented a fifth order Bessel lowpass filter. 
The poles of a normalized filter (o) 3 dB= 1 rad/s) are located at si, 2 = 
-0.9576766 ±j1 .471 1244, s 3 , 4 =-1.3808774 ±j0.7179096, and S5 = 

-1.5023160 [9]. The normalized filter system function can be 
written as 


H 0 (s) = = _dl_ + _d2 ^L_ 

(s- Ji)(.y- s 2 )•* •(£-•%) s-s 1 s-s 2 s-s s 


(17) 


The system function of the actual filter can be obtained by the 
frequency transformation H(s)=H 0 (s/co c ) with cd c the 3 dB cutoff 
frequency. The impulse response is obtained by taking the inverse 
Laplace transform of H(s), as 

h c (t) = \co c e w +A 2 ( 0 c e S2t0c, +-"+A 5 0 ) c e s>,Oct (18) 

The values of Aj's can be obtained through partial fraction, as 
Ai ,2 = 0.6193736 ±j1. 0986409, A 3>4 =-4.9168825 ±j(-2.61 31 445), 
and A 5 = 8.5950178. 


5. The Computer Program 

Figure 2 shows a flowchart of the program. It is written in 
Fortran 77. The Poisson and Gaussian random number generators are 
from [4], The uniform random variable generator is also from [4], 
Other subroutines are developed by ourselves. 

Figures 3 shows sample histograms of the Gaussian, the Poisson, the 
McIntyre, and the Webb random variable generators. They all agreed 
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very well with the corresponding theoretical distributions. The 
fluctuations at the tails of the histograms were due to the finite 
sample sizes. 

Figure 4 shows the impulses of two low pass filters, a fifth order 
Bessel lowpass filter with 3 dB bandwidth equal to 100 MHz and an 
ideal integrator with integration time equal to 10 ns. Figures 5 and 
6 show plots of the sample output from this APD preamplifier 
simulator using the lowpass filter given in Figure 4. The input 
optical signal had a rectangular pulse shape of width 10 ns. The 
received background radiation power was 1.0 nW and the peak signal 
power including the background was 20 nW. The pulse energy was 
therefore 1067 photons, or 373 primary photoelectron per pulse 
assuming 35% quantum efficiency. Other parameter values used are 
listed in Table 1. 
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Figure 1 . Block diagram of the APD preamplifier simulator 
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Figure 2. Block diagram of the APD preamplifier simulator program. The 
names of the subroutines and sub-subroutines are given in parenthesis. 
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Histogram of the Gaussian R.V. Generator 
for jt-0 and <M.O 



Output 


Histogram of the Poisson R.V. Generator 
for nbar-5 
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Histogram of the McIntyre A.V. Generator 
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Number of Secondary Photoeiectron* 


Histogram of the Webb's R.V. Generator 

«v*r«o« primary photoalactrons-10, G-200, and sample ***-100,000 



Number of Secondary Photoelectrons 


Figure 3. Sample histograms of the r.v. generators. 
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Impulse Response of the Lowpass Filter 



time (ns) 


Figure 4. Impulse responses of a fifth order Bessel lowpass filter 
and an ideal integrator of integration time Ax=10 ns. 
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APD Preamp Simulator Output 



t (ns) 


Figure 5. Sample APD preamp and lowpass filter outputs. The filter 
used was a fifth order Bessel lowpass filter with f3dB=100 MHz. 
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APD Preamp Simulator Output 
Integrator ( x=1 0 ns) as LPF 



t (ns) 


Figure 6. Sample APD preamp and lowpass filter outputs. The filter 
used was an ideal integrator with integration time Ax=10 ns. 
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Table 1. System Parameter Values for Figure 4. 


laser pulse shape rectangular 

laser pulsewidth 10 ns 

laser pulse peak power 1 9 nW 

background radiation power 1.0 nW 

wavelength, lambda 1060 nm 

APD quantum efficiency, eta 35% 

average APD gain, G 2 0 0 

APD ionization coefficient ratio, k e ff 0.010 

APD surface leakage current, Is 15 nA 

APD bulk leakage current, lb 50 pA 

preamplifier feedback resistor, R 20 kohm 

preamplifier noise temperature, Tn 700 0 K 


noise filter type 5th order Bessel LPF 

noise filter 3dB bandwidth, fc 100 MHz 


# points sampled, N 
sampling interval, dt 


2048 
0.10 ns 



makefile 
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apdout: apdouto POIDEVax) GAMMLNa.o RANI ax> mddev.o webdev.o GASDEVa.o \ 

fltresp.o Ipfout.o 

f77 -o apdout apdouto POIDEVa.o GAMMLNa.o RAN1a.o mddev .0 webdev.o \ 
GASDEVa.o fltresp.o Ipfout.o 
POIDEVao: POIDEVaf 

f77 -c POIDEVa.f 
GAMMLNao: GAMMLNa.f 

f77 -c GAMMLNa.f 
RANIao: RANla.f 

f77 -c RANla.f 
mcidev.o: mcidev.f 

f77 -c mcidev.f 
webdev.o: webdev.f 

f77 -c webdev.f 
GASDEVao: webdev.f 

f77 -c GASDEVa.f 
fltresp.o: fltresp.f 

f77 -c fltresp.f 
Ipfout.o: Ipfout.f 

f77 -c Ipfout.f 
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apdout.f 


* This program generate an array of random variable to simulate 

* the APD and the preamp output equivalent at the input to the 

* preamp (in ampere) in response to an input array of samples of 

* incident optical power in watts 


c 

INTEGER N, I, idum, m, II, nl, kk, NLPF 
PARAMETER (N-40000.NLPF-1000) 

REAL lambda, eta, G, Keff, Is, lb, R, Tn, q 
PARAMETER (q-1.6e-19) 

REAL dt, nbar 
REAL hf, F, sigma, mu, xl 
REAL PS(N), Y(N), Y1(N), RESPNS(NLPF) 
c 

c Received optical signal and background radiation power 
c 

DATA PS / 500*1.0e-9,1 00*2. Oe-8, 500*1. 06-9,8900*0.0, 30000*3.36-9 / 
c 

REAL mcidev, poideva, gasdeva, webdev (subroutines to be called 
c 

c System Parameters and definitions of some of the variables used 
c 

c lambda-laser diode wavelength (nm) 
c eta-APD quantum efficiency 
c G-average APD gain 

c Keff-ratio of ionization coefficients of APD 
c Is-APD surface leakage current in nanoamperes (A) 
c Ib-APD bulk leakage current in nanoamperes (A) 
c R- Preamp transimpedance (APD load resistance) (Ohm) 
c Tn-equivalent noise temperature (Kalvin) 
c q-electron charge (C) 
c dt-sampling interval (s) 
c nbar-average # of primary photoelectrons in dt 
c PS(l)=input optical signal array 
c N-number of elements in PS(I) 
c 

lambda-1060.0 
eta-0.35 
G-200 
Keff-0.010 
ls=15.0e-9 
lb-0.050e-9 
R-20000.0 
Tn-700.0 
dt-0.1 0e-9 
c 

print *, 'Seed for the random number generator?' 
read *, idum 
c 

c Other parameters used in the calculation 
c 

hf-1 242.0/lambda*q I photon energy 
F-Keff*G+(2.0-1.0/G)*(1.0-Keff) I APD excess noise factor 
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sigma-sqrt(2.0*1.38e-23*Tn/R/dt+q*ls/dt) ! rms Gaussian noise 
mu-ls ! mean of the Gaussian noise 

c 

c Start calculation 
c 

do 500 1-1 ,N 
c 

c Number of secondary photoelectrons from the APD 
c 

nbar-eta*PS(l)/hf*dt+lb*dt/q 

m-0 

if( nbar.LE.10.0 ) then 

nl- poideva(nbar.idum) 
if(nl.EQ.O) go to 200 
do 100 11-1, nl 

100 m-m+mcidev(G,Keff,idum) 

else 

m-webdev(nbar,G,F,idum) 

end if 
c 

c Preamp Gaussian noise 
c 

200 x1=sigma*gasdeva(idum)+mu 
c 

c Total output voltage (volts) 
c 

Y(l)=x1+real(m)*q/dt ! current at the input to the preamp 
Y(I)=Y(I)*R I output voltage from preamp 

500 continue 
c 

c The digital lowpass filter 
c 

fc=1 00.0e6 I the 3dB cutoff frequency (Hz) 

call fltresp(NLPF, fc, dt, RESPNS) ! Lowpass filter impulse response 
call lpfout(Y, N, RESPNS, NLPF, Y1) I convolution 
c 

c Print the Results 
c 

write(20,700) 

700 format(1x,'index’,7x,'t(ns)',7x,'lnput(nW) , ,5x, 

+ ’APDoutputtVy^x.’LPFoutpuUV)' ) 

do 900 kk-1 ,N 

write(20,800) kk, dt*kk*1.0e9, PS(kk)*1 .0e9, Y(kk), Y1(kk) 
800 format(1x, 16, 2f14.3, 2e14.3 ) 

900 continue 
c 

stop 

end 
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FUNCTION poideva(xm,idum) 

INTEGER idum 
REAL poideva.xm.PI 
PARAMETER (Pl=3. 141 592654) 
c 

REAL alxm,em,g,oldm,sq,t,y ,gammlna,ran 1 a 
SAVE alxm,g,oldm,sq 
DATA oldm /-I./ 
if (xm.lt.12.)then 
if (xm.ne.oldm) then 
oldm=xm 
g=exp(-xm) 
endif 
em=-1 
t=1. 

2 em=em+1. 

t=t*ran1a(idum) 
if (t.gt.g) goto 2 
else 

if (xm.ne.oldm) then 
oldm=xm 
sq=sqrt(2.*xm) 
alxm=log(xm) 

g=xm*alxm-gammlna(xm+1 .) 
endif 

1 y=tan(PI*ran1a(idum)) 

em=sq*y+xm 
if (em.lt.O.) goto 1 
em=int(em) 

t=0.9*(1 .+y**2)*exp(em*alxm-gammlna(em+1 .)-g) 
if (ranla(idum).gt.t) goto 1 
endif 

poideva=em 

return 

END 
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GAMMLNa.f 


FUNCTION gammlna(xx) 

REAL gammlna.xx 
INTEGER j 

DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) 

SAVE cof.stp 

DATA cof,stp/76. 1 80091 729471 46d0, -86.50532032941 677d0, 
*24.01409824083091 d0,- 1 .231 739572450 155d0,. 12086509738661 79d-2, 
*-. 5395239384953d-5, 2. 506628274631 0005d0/ 

x=xx 

y=x 

tmp=x+5.5d0 

tmp=(x+0.5d0)*log(tmp)-tmp 
ser=1 .00000000019001 5d0 
do 11 j=1,6 
y=y+1.d0 
ser=ser+cof(j)/y 
1 1 continue 

gammlna=tmp+log(stp*ser/x) 

return 

END 
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RANla.f 


FUNCTION ranla(idum) 

INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV 
REAL ranla.AM.EPS.RNMX 

PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, 
*NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1 .2e-7,RNMX=1 .-EPS) 

INTEGER j,k,iv(NTAB),iy 
SAVE iv.iy 

DATA iv /NTAB*0/, iy 101 
if (idum.le.O.or.iy.eq.O) then 
idum=max(-idum,1) 
do 11 j=NTAB+8,1,-1 
k=idum/IQ 

idum=IA‘ r (idum-k*IQ)-l R*k 
if (idum.lt.O) idum=idum+IM 
if (j.le.NTAB) iv(j)=idum 
1 1 continue 
iy=iv(1) 
endif 

k=idum/IQ 

idum=IA*(idum-k*IQ)-IR*k 
if (idum.lt.O) idum=idum+IM 
j=1+iy/NDIV 
jy=iv(j) 
iv(j)=idum 

ran1a=min(AM*iy,RNMX) 

return 

END 
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mcidev.f 


This subroutine generates a random number according to the * 

Mclntyre-Conradi distribution for the number of * 

secondary electrons output from an APD in response to one * 

* primary photoelectron. 

*************’**'*** , *'******'*'*****'*'* , 0 #* , ****'**»*##'*'********.*^^^^.^^ 

c 

FUNCTION mcidev(G,k,idum) 

INTEGER idum 
REAL mcidev,G,k,PI 
PARAMETER (Pl=3. 141 592654) 
c 

REAL aO,xO,cO,y,t 
REAL ranla.pmcint.cmcint 
a0=G/3.0 
x0=-G/50.0 
c0=4.0/G 

1 y=aO*tan(PI*ran1a(idum))+xO 

if (y.lt.O.) goto 1 
y=int(y) 

t=pmcint(y,G,k)/cmcint(y,aO,xO,cO) 
if (ranla(idum).gt.t) goto 1 
mcidev=y 
return 
END 
c 
c 

c The McIntyre distribution 
c 

FUNCTION pmcint(y,G,k) 

REAL pmcint.y.G.k 
REAL k1,G1 

DOUBLE PRECISION x1,x2,x3,x4 
INTEGER m,i 
c 

kl = 1 ,0-k 
G1=G-1 .0 
m = int(y-2.0) 
c 

x 1 =log( (1+k*G1)/G ) 
x2 = 1.0 + k*y/k1 
x3=log( G1/G ) 
x4=0.0 
do 100 i=0, m 

100 x4= x4 + log( 1.0+ k*real(i)/( y-real(i) ) ) 

pmcint= exp( x1*x2+x3*(y-1.0)+x4 ) 
return 
end 
c 

c The comparison function (Lorentzian) 
c 

FUNCTION cmcint(x,a0,x0,c0) 

REAL cmcint.x.aO.xO.cO 

cmcint=aO*aO*cO/( a0*a0+(x-x0)*(x-x0) ) 

return 
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webdev.f 


This subroutine generates a random number according to the * 
* Webb's distribution for the number of * 

secondary electrons output from an APD in response to n=10 * 

primary photoelectrons. * 

c 

FUNCTION webdev(nbar,G,F,idum) 

INTEGER idum 
REAL webdev,nbar,G,F,PI 
PARAMETER (Pl=3. 141 592654) 
c 

REAL aO,xO,cO,y,t 
REAL ranla.pwebb.cwebb 

a0=G*sqrt(2.0*F*nbar) 

xO=G*nbar 

c0=1 .0/(G*sqrt(F*nbar)) 

1 y=aO*tan(PI*ran1a(idum))+xO 

if (y.lt.O.) goto 1 
y=int(y) 

t=pwebb(y,nbar,G,F)/cwebb(y,aO,xO,cO) 
if (ranla(idum).gt.t) goto 1 
webdev=y 
return 
END 
c 
c 

c The Webb distribution 
c 

FUNCTION pwebb(y,nbar,G,F) 

REAL pwebb,y,G,F,nbar,PI 
PARAMETER (Pl=3. 141592654) 

REAL x1,x2,x3,x4 
c 

x1=2.0*G*G‘F*nbar 

x2=y-G*nbar 

x3=1.0 + (y-G*nbar)*(F-1.0)/(G*F*nbar) 
x4=1.0/sqrt(PI*x1)/x3/sqrt(x3) 
pwebb=x4*exp(-x2*x2/x1/x3) 
return 
end 
c 

c The comparison function (Lorentzian) 
c 

FUNCTION cwebb(x,aO,xO,cO) 

REAL x,aO,xO,cO,cwebb 

cwebb=aO*aO*cO/( aO*aO+(x-xO)*(x-xO) ) 

return 

end 
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GASDEVa.f 


FUNCTION gasdeva(idum) 
INTEGER idum 
REAL gasdeva 
CU USES rani a 
INTEGER iset 

REAL fac,gset,rsq,v1 1 v2,ran1a 
SAVE iset,gset 
DATA iset/0/ 
if (iset.eq.O) then 
1 v1=2.*ran1a(idum)-1. 

v2=2.*ran1a(idum)-1. 
rsq=v1 **2+v2**2 
if(rsq.ge.1..or.rsq.eq.0.)goto 1 
fac=sqrt(-2.*log(rsq)/rsq) 
gset=v1 *fac 
gasdeva=v2*fac 
iset=1 
else 

gasdeva=gset 

iset=0 

endif 

return 

END 
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fltresp.f 


subroutine fltresp(M,fc,Delt,h) 
integer M, k 
real h(M) 
real fc, Delt, he 
c 

do 100 k=1,M 
h(k)=Delt*hc(k*Delt,fc) 

100 continue 
return 
end 
c 

c The continous time impulse response function 
c 

function hc(t,fc) 
real hc,t,fc 
integer M1,k 

parameter (M1=5,PI=3.1416) 
complex A(M1), s(M1), chc 
c 

A(1) = (0. 6 163736, 1.0986409) 

A(2) = (0.61 63736,-1 .0986409) 
A(3) = (-4.91 68825,-2.61 3 1445) 
A(4) = (-4.91 68825,2.61 31 445) 
A(5) = (8. 5950178, 0.0) 
c 

s(1) = (-0. 9576766,1 .471 1244) 
s(2) = (-0.9576766, -1.471 1244) 
s(3) = (-1 .3808744,0.7179096) 

S(4) = (-1 .3808744,-0.7179096) 
s(5) = (-1 .5023160,0.0) 
c 

hc=0.0 

do 200 k=1,M1 

chc=A(k)*cexp( s(k)*2.0*PI*fc*t ) 
200 hc=hc+real(chc) 
hc=hc*2.0*PI*fc 
return 
end 
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Ipfout.f 


A subroutine to perform convolution of array 'data’ with * 
array 'respns'. The dimension of the first array, nd, is 
assumed to be greater than the second, nresp. The resultant * 
array is given by 'output*. Zero padding was used when * 
calculating the first nresp elements. 


subroutine lpfout(data,nd,respns,nresp,output) 

integer nd, nresp, i, j, jO, jj 

real data(nd), respns(nresp), output(nd) 

do 100 i-1,nd 
output(i)-0.0 

if(i.GE.nresp) then 
j0-i-nresp+1 
do 20 H0,i 
jj-i-j + 1 

20 output(i)=output{i)+data(j)*respns(jj) 
else 

do 40 j=1,i 
jj-i-j + 1 

4 0 output (i)=output(i)+data(j)*respns(jj) 

end if 

1 00 continue 
return 
end 



